Oscillons¶

In [1]:
%matplotlib inline
In [2]:
import os,sys
sys.path.append('./bubbles_codes/')
from plotting import *
from bubble_tools import *
from experiment import *
%run './bubbles_codes/plotting.py'
from itertools import chain
import operator as op
import pywt
In [3]:
def get_bubble(exp_params, sim, crit_thresh, crit_rad):
    path2RESTsim  = rest_sim_location(*exp_params, sim)
    path2CLEANsim = clean_sim_location(*exp_params, sim)

    bubble = np.load(path2CLEANsim)[0]
    bubble = np.array([bubble[0]])
    totbeta = np.load(path2RESTsim)[-1]

  #  bubble = multiply_bubble(bubble, lightc, phieq, totbeta, normal, nLat)
    bubble = bubble[0]
  #  bubble = gaussian_filter(bubble, sigma=2, mode='nearest')

    tcen, xcen = find_nucleation_center(bubble, phieq, crit_thresh, crit_rad)
    bubble = bubble[max(tcen-nLat*2//3,0):tcen]
    mn = np.mean(bubble)
    bubble-= mn
    return bubble, mn

def get_HT(array):  
    array = array - np.mean(array)
    w  = np.fft.fftfreq(len(array), d=1) * len(array)
    FD = np.fft.fft(array, axis=0)
    FD[w<=0.,:] = 0.
    FD[w>0.,:] *= 2.
    HD = np.fft.ifft(FD, axis=0)
    return np.abs(HD)

def get_4T(array):
    nT, nN = np.shape(array)
    tw = np.fft.fftfreq(nT, d=1) * nT
    xw = np.fft.fftfreq(nN, d=1) * nN

    tFD = np.fft.fft(array, axis=0)
    tFD[tw<0,:] = 0.
    tFD[tw>0,:]*= 2.
    tHD = np.fft.ifft(tFD, axis=0)

    xFD = np.fft.fft(array, axis=1)
    xFD[:,xw<0] = 0.
    xFD[:,xw>0]*= 2.
    xHD = np.fft.ifft(xFD, axis=1)

    lFD = np.fft.fft(tHD, axis=1)
    lFD[:,xw<0] = 0.
    lFD[:,xw>0]*= 2.
    lHD = np.fft.ifft(lFD, axis=1)

    rFD = np.fft.fft(xHD, axis=0)
    rFD[tw>0,:] = 0.
    rFD[tw<0,:]*= 2.
    rHD = np.fft.ifft(rFD, axis=0)
    return np.abs(tHD), np.abs(xHD), np.abs(lHD), np.abs(rHD)

def get_cummul_power(freqsT, freqsX, PST, PSX):
    # flat list of all frequencies where response is strong
    all_frqT = flatten_comprehension2(savefrqT)
    all_frqX = flatten_comprehension2(savefrqX)

    all_psT  = flatten_comprehension(savePST)
    all_psX  = flatten_comprehension(savePSX)

    # cut duplicates
    set_Tfreqs = np.sort(all_frqT)
    set_Xfreqs = np.sort(all_frqX)

    # add field responses to get relative importance of each mode
    set_Ttotamps = np.zeros(len(set_Tfreqs))
    for fi, frq in enumerate(set_Tfreqs):
        inds = np.argwhere(np.abs(all_frqT)==frq).flatten()
        set_Ttotamps[fi] = np.sum(all_psT[inds])

    # and in X
    set_Xtotamps = np.zeros(len(set_Xfreqs))
    for fi, frq in enumerate(set_Xfreqs):
        inds = np.argwhere(np.abs(all_frqX)==frq).flatten()
        set_Xtotamps[fi] = np.sum(all_psX[inds])
    return set_Tfreqs, set_Xfreqs, set_Ttotamps, set_Xtotamps

def get_osc_trajectory(array, extent):
    # find maximum value in array
    T, N = np.shape(array)
    flattened = chain.from_iterable(array)
    max_idx, max_val = max(enumerate(flattened), key=op.itemgetter(1))
    row = max_idx // N
    col = max_idx % N
#    print(array[row][col] == max_val, max_val, array[row][col], row, col)
#    ax[axis].plot(col, row, 'ro', ms=10)

    # use as starting point for oscillon trajectory
#    print(col, row)
    maxLine = []
    CCol = col
    for rr in range(row)[::-1]:
        if rr == row-1:
            col = CCol
        colmin = col - extent
        colmax = col + extent+1
        val = 0
        for cc in range(colmin, colmax):
            cc = cc%N
            if array[rr][cc] > val:
                val = array[rr][cc]
                col = cc
        maxLine.append(col)
    maxLine = maxLine[::-1]
    for rr in range(row, T):
        if rr == row:
            col = CCol
        colmin = col - extent
        colmax = col + extent+1
        val = 0
        for cc in range(colmin, colmax):
            cc = cc%N
            if array[rr][cc] > val:
                val = array[rr][cc]
                col = cc
        maxLine.append(col)
    return np.array(maxLine)

# average oscillon trajectories
def tolerant_mean(arrs):
    lens = np.array([len(i) for i in arrs])
  #  print(lens.tolist())
    arr  = np.zeros((len(lens), np.max(lens)))
    for ri, osc in enumerate(arrs):
        arr[ri, :len(osc)] = osc
    arr[arr == 0] = np.nan
    return np.nanmean(arr, axis=0), np.nanstd(arr, axis=0)

def tolerant_mean2d(arrs):
    lens = np.array([np.shape(i) for i in arrs])
    lens1, lens2 = lens[:,0], lens[:,1]

    arr  = np.zeros((len(arrs), np.max(lens1), np.max(lens2)))
    for ri, osc in enumerate(arrs):
        n1, n2 = np.shape(osc)
        arr[ri, :n1, :n2] = osc
    arr[arr == 0.] = np.nan
    return np.nanmean(arr, axis=0), np.nanstd(arr, axis=0)

def flatten_comprehension(matrix):
    return np.array([item for row in matrix for item in row])

def flatten_comprehension2(matrix):
    return np.array([np.abs(np.round(item,10)) for row in matrix for item in row])
In [4]:
simLists = []
tmp = 0
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
for sim in range(minSim, maxSim):
    path2RESTsim = rest_sim_location(*exp_params, sim)
    if os.path.exists(path2RESTsim):
        path2CLEANsim = clean_sim_location(*exp_params, sim)
        if os.path.exists(path2CLEANsim):
            simLists.append(sim)
In [5]:
crit_rad = 10
crit_thresh = right_Vmax.x + 4.*sigmafld
print(crit_thresh)

deltae = 25
print(deltae * dx * np.sqrt(m2(lamb)))
5.348981768243531
1.9301011109426147
In [6]:
saveData = False
In [7]:
if saveData:
    oscillonData = []
    for si, sim in enumerate(simLists):
        bubble, mn = get_bubble(exp_params, sim, crit_thresh, crit_rad)
        nT, nN = np.shape(bubble)

        try: aHT = get_HT(bubble)
        except: continue

        maxln1 = get_osc_trajectory(aHT, extent=deltae)
        avOscillon = bubble[np.arange(nT), maxln1][::-1]
        maxln1-= nN//2

        oscillonData.append(np.array([sim, avOscillon, mn]))

    np.save('./oscillon_data.npy', oscillonData)

else:
    oscillonData = np.load('./oscillon_data.npy')
    print(np.shape(oscillonData))
    simLists, oscillonAmp, means = oscillonData[:,0], oscillonData[:,1], oscillonData[:,2]

    oscillonFreqData = []
    for si, sim in enumerate(simLists):
        oscillon = oscillonAmp[si]
        mn = means[si]

        nT = len(oscillon)
        fft_oscillon = np.fft.rfft(oscillon)
        freqs = np.fft.rfftfreq(nT, d=1) * nT
        PS_oscillon = (np.conj(fft_oscillon)*fft_oscillon).real
        ind = np.argmax(PS_oscillon)

        oscillonFreqData.append(np.array([freqs, PS_oscillon, ind]))
    oscillonFreqData = np.array(oscillonFreqData)
(653, 3)
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:34: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
In [8]:
# One example
sim = simLists[0]
bubble, _ = get_bubble(exp_params, sim, crit_thresh, crit_rad)
nT, nN = np.shape(bubble); print(nT, nN)

aHT = get_HT(bubble)
tHD, xHD, lHD, rHD = get_4T(bubble)

fig, ax = plt.subplots(3, 2, figsize=(6,6))
im0 = ax[0,0].imshow(bubble, aspect='auto', interpolation='none', origin='lower')
im1 = ax[0,1].imshow(aHT, aspect='auto', interpolation='none', origin='lower')
im2 = ax[1,0].imshow(tHD, aspect='auto', interpolation='none', origin='lower')
im3 = ax[1,1].imshow(xHD, aspect='auto', interpolation='none', origin='lower')
im4 = ax[2,0].imshow(lHD, aspect='auto', interpolation='none', origin='lower')
im5 = ax[2,1].imshow(rHD, aspect='auto', interpolation='none', origin='lower')

im = [im0,im1,im2,im3,im4,im5]
titls = [r'$\rm Original$', r'$\rm Hilbert \; T$', r'$\rm Hilbert \; T$', 
         r'$\rm Hilbert \; X$', r'$\rm Left \; travelling$', r'$\rm Right \; travelling$']
for ii in range(len(im)):
    ax[np.divmod(ii,2)].set_title(titls[ii])
    plt.colorbar(im[ii], ax=ax[np.divmod(ii,2)])
plt.tight_layout()
plt.show()
682 1024
In [9]:
nT, nN = np.shape(bubble)
aHD = get_HT(bubble)

amp = 3
t = aHD**amp / np.sum(aHD**amp, axis=1)[:,None] 
SB = t*np.arange(nN)[None,:]
xt = np.sum(SB, axis=1)

fig, ax = plt.subplots(1, 3, figsize=(12,3))
im1 = ax[0].imshow(bubble, aspect='auto', interpolation='none', origin='lower') #original
im2 = ax[1].imshow(aHD,     aspect='auto', interpolation='none', origin='lower') #hilbert transform amplitude
im3 = ax[2].imshow(SB,      aspect='auto', interpolation='none', origin='lower') #power spectrum method

l1 = get_osc_trajectory(aHD, extent=deltae)
l2 = gaussian_filter1d(l1, sigma=10, mode='nearest')

l3 = get_osc_trajectory(SB, extent=deltae)
l4 = gaussian_filter1d(l3, sigma=10, mode='nearest')

ylist = np.arange(nT)
ax[1].plot(l1, ylist, color='m', linewidth=3)
ax[1].plot(l2, ylist, color='orange', linewidth=1)

ax[2].plot(l3, ylist, color='white', linewidth=3)
ax[2].plot(l4, ylist, color='orange', linewidth=1)

plt.colorbar(im1, ax = ax[0])
plt.colorbar(im2, ax = ax[1])
plt.colorbar(im3, ax = ax[2])
plt.tight_layout()
plt.show()
In [10]:
# wavelet transform on oscillon maplitude
C,B    = 1,1 # center, bandwidth # C = freq to investigate, B = bandwidth
wavelt = "cmor%f_%f"%(C,B)
dt     = 1
scales = np.arange(1, nT+1)

# construct field value of oscillon over trajectory above
oscillonAmp = bubble[np.arange(nT), l1] 
CC, f = pywt.cwt(oscillonAmp, scales, wavelt, sampling_period=dt)

smoscillonAmp = bubble[np.arange(nT), l2]
BB, f = pywt.cwt(smoscillonAmp, scales, wavelt, sampling_period=dt)

fig, ax = plt.subplots(3, 1, figsize=(5, 5), gridspec_kw=dict(hspace=0), sharex=True)
ax[0].imshow(np.abs(CC), aspect='auto', interpolation='none', origin='lower')
ax[1].imshow(np.abs(BB), aspect='auto', interpolation='none', origin='lower')
ax[2].plot(oscillonAmp)
ax[2].plot(smoscillonAmp)
plt.tight_layout()
plt.show()
In [11]:
#oscillonWT = []
if not saveData:
    for si, sim in enumerate(simLists):
        if si%10!=0: continue

        bubble, _ = get_bubble(exp_params, sim, crit_thresh, crit_rad)
        nT, nN = np.shape(bubble)
        try: aHD = get_HT(bubble)
        except: continue

        maxln1 = get_osc_trajectory(aHD, extent=deltae)
        avOscillon = bubble[np.arange(nT), maxln1][::-1]
        maxln1-= nN//2

    #    # wavelet transform on oscillon maplitude
    #    C,B     = 1,1 # center, bandwidth # C = freq to investigate, B = bandwidth
    #    wavelet = "cmor%f_%f"%(C,B)
    #    dt      = 1
    #    scales  = np.arange(1,nT)
    #    CC, f   = pywt.cwt(avOscillon, scales, wavelet, sampling_period=dt)
    #    oscillonWT.append(CC)

        print('SIMULATION', sim)
        fig, ax = plt.subplots(1, 3, figsize=(9, 3))

        exts = np.array([-nN//2,nN//2,0,nT]) * dx * np.sqrt(m2(lamb))
        im1 = ax[0].imshow(bubble, extent=exts, aspect='auto', interpolation='none', origin='lower', cmap='RdBu') #original

        ylist = np.arange(nT) * dx * np.sqrt(m2(lamb))
        ax[0].plot(maxln1 * dx * np.sqrt(m2(lamb)), ylist, ls='-', color='k', linewidth=1)

        ax[0].set(xlabel=r'$m \, r$')
        ax[0].set(ylabel=r'$m \, t$')
        ax[0].tick_params(direction='in', which='both', top=True, right=True)
        ax[0].grid(ls=':', color='darkgray', alpha=0.5)

        cbar = fig.colorbar(im1, ax=ax[0], ticks=mticker.MultipleLocator(np.pi/2), \
                            format=mticker.FuncFormatter(multiple_formatter()))
        cbar.ax.set_title(r'$\bar{\phi}$')

        ax[1].plot(np.arange(nT), avOscillon)
        ax[1].set_xlabel(r'$t$')
        ax[1].set_ylabel(r'$\bar{\phi}(osc,t)$')
        ax[1].axhline(0, ls=':', color='k')

        y_fft = np.fft.rfft(avOscillon)
        freqs = np.fft.rfftfreq(nT, d=1) * nT
        psyfd = np.conj(y_fft)*y_fft
        ind   = np.argmax(psyfd)

        ax[2].loglog(freqs[ind], psyfd[ind], 'o')
        ax[2].loglog(freqs, psyfd, linewidth=0.5)
        ax[2].set_xlabel(r'$k$')
        ax[2].set_ylabel(r'$P(k)$')

        for aa in ax: aa.grid(True)

        plt.tight_layout()
        plt.show()
SIMULATION 7
/home/dpirvu/.local/lib/python3.7/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part
  return np.asarray(x, float)
SIMULATION 60
SIMULATION 122
SIMULATION 185
SIMULATION 255
SIMULATION 316
SIMULATION 369
SIMULATION 421
SIMULATION 475
SIMULATION 521
SIMULATION 583
SIMULATION 608
SIMULATION 651
SIMULATION 685
SIMULATION 769
SIMULATION 835
SIMULATION 901
SIMULATION 974
SIMULATION 1036
SIMULATION 1097
SIMULATION 1164
SIMULATION 1249
SIMULATION 1316
SIMULATION 1432
SIMULATION 1493
SIMULATION 1539
SIMULATION 1589
SIMULATION 1625
SIMULATION 1698
SIMULATION 1755
SIMULATION 1845
SIMULATION 1942
SIMULATION 2023
SIMULATION 2081
SIMULATION 2111
SIMULATION 2201
SIMULATION 2252
SIMULATION 2285
SIMULATION 2331
SIMULATION 2402
SIMULATION 2468
SIMULATION 2533
SIMULATION 2573
SIMULATION 2612
SIMULATION 2724
SIMULATION 2775
SIMULATION 2824
SIMULATION 2903
SIMULATION 2961
SIMULATION 3019
SIMULATION 3118
SIMULATION 3155
SIMULATION 3233
SIMULATION 3305
SIMULATION 3370
SIMULATION 3432
SIMULATION 3481
SIMULATION 3528
SIMULATION 3607
SIMULATION 3646
SIMULATION 3696
SIMULATION 3759
SIMULATION 3800
SIMULATION 3866
SIMULATION 3926
SIMULATION 3987

FT of Oscillon Signal¶

In [12]:
if saveData:

    savefrqT, savefrqX, savePST, savePSX = [], [], [], []
    for si, sim in enumerate(simLists):
        bubble, _ = get_bubble(exp_params, sim, crit_thresh, crit_rad)
        nT, nN = np.shape(bubble)

        bubbleT = np.fft.fft(bubble, axis=0)
        freqsT  = np.fft.fftfreq(nT, d=1) * nT # get frequencies; d = x space units
        freqsT  = np.outer(freqsT, np.ones(nN))
        PST     = (np.conj(bubbleT)*bubbleT).real

        bubbleX = np.fft.fft(bubble, axis=1)
        freqsX  = np.fft.fftfreq(nN, d=1) * nN
        freqsX  = np.outer(np.ones(nT), freqsX)
        PSX     = (np.conj(bubbleX)*bubbleX).real

        valMax  = 2e4
        coordsT = (PST>=valMax) # separate freqs where there is large field response
        coordsX = (PSX>=valMax)

        savefrqT.append(freqsT[coordsT])
        savefrqX.append(freqsX[coordsX])

        savePST.append(PST[coordsT])
        savePSX.append(PSX[coordsX])

    set_Tfreqs, set_Xfreqs, set_Ttotamps, set_Xtotamps = get_cummul_power(freqsT, freqsX, PST, PSX)
    del freqsT, freqsX, PST, PSX

    np.save('./oscillon_freqs.npy', np.array([set_Tfreqs, set_Xfreqs, set_Ttotamps, set_Xtotamps]))
In [13]:
# plot frequency vs response; plot mass from frequency and compare with free mass
set_Tfreqs, set_Xfreqs, set_Ttotamps, set_Xtotamps = np.load('./oscillon_freqs.npy')
masses = set_Tfreqs[:,None]**2. - set_Xfreqs[None,:]**2.
masses = masses.flatten()
masses = masses[masses>0.]

fig, ax = plt.subplots(1, 3, figsize=(9, 3))
ax[0].plot(set_Tfreqs, set_Ttotamps)
ax[0].set_xlabel(r'$\omega$')
ax[1].plot(set_Xfreqs, set_Xtotamps)
ax[1].set_xlabel(r'$k$')
ah = ax[2].hist(masses, bins=10, density=True)
ax[2].set_xlabel(r'$m$')

#for aa in ax: aa.set_yscale('log')
plt.tight_layout()
plt.show()
In [14]:
def omega(ksq, a, m):
    return a*ksq + m**2

if False:
    fig, ax = plt.subplots(1, 2, figsize=(6, 2))
    for si, sim in enumerate(simLists):
        bubble, _ = get_bubble(exp_params, sim, crit_thresh, crit_rad)
        nT, nN = np.shape(bubble)

        bubblePS = np.fft.fft2(bubble)[:nT//6,:nN//6]
        bubblePS = np.conj(bubblePS)*bubblePS

        maxVal = 5e6
        bubble_coords = np.argwhere(bubblePS>maxVal)

        kkk, ooo = bubble_coords[:,1]**2., bubble_coords[:,0]**2.
        popt2, pcov2 = sco.curve_fit(omega, kkk, ooo)

        if True:
            if si%10!=4: continue
            ax[0].imshow(bubblePS.real>maxVal, aspect='auto', interpolation='none', origin='lower')

            ax[1].plot(kkk, omega(kkk, *popt2), '-')#, label='fit: a=%5.3f, m=%5.3f' % tuple(popt2))
            ax[1].plot(kkk, ooo, 'go', ms=1)
    for aa in ax:
        aa.set_xlabel(r'$k^2$')
        aa.set_ylabel(r'$\omega_k^2$')
    plt.show()

Average Oscillon Signal¶

In [15]:
fig, ax = plt.subplots(1, 2, figsize = (8, 2.5))

simLists, oscillonAmp, means = oscillonData[:,0], oscillonData[:,1], oscillonData[:,2]
for si, sim in enumerate(simLists):
    oscillon = oscillonAmp[si]
    mn = means[si]
    freqs, PS_oscillon, ind = oscillonFreqData[si]

    nT = len(oscillon)
    ax[0].plot(np.sqrt(m2(lamb))*dx*np.arange(nT), oscillon + mn, linewidth=0.05)

    ax[1].loglog(freqs[ind], PS_oscillon[ind], 'o', ms=1)
    ax[1].loglog(freqs, PS_oscillon, linewidth=0.05)

avOS, error1 = tolerant_mean(oscillonAmp)
ax[0].plot(np.sqrt(m2(lamb))*dx*np.arange(len(avOS)), avOS + np.mean(means), linewidth=1, color='k')

avSP, error1 = tolerant_mean(oscillonFreqData[:,1])
avfq, error1 = tolerant_mean(oscillonFreqData[:,0])
ax[1].plot(avfq, avSP, linewidth=1, color='k')

ax[0].yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))
ax[0].yaxis.set_major_formatter(plt.FuncFormatter(multiple_formatter()))
ax[0].axhline(phieq, ls=':', color='k', linewidth=1)

ax[0].set_xlabel(r'$m \, t$')
ax[0].set_ylabel(r'$\bar{\phi}^{\rm osc}(t)$')

ax[0].set_xlim(0, 18)#np.sqrt(m2(lamb))*dx*len(avOS))

ax[1].set_xlabel(r'$k$')
ax[1].set_ylabel(r'$P(k)$')

for aa in ax: aa.grid(True, ls=':', alpha=0.5)
plt.tight_layout()
plt.savefig('./plots/av_oscillon_and_freq.pdf', dpi=500, rasterize=True)
plt.show()
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:128: RuntimeWarning: Mean of empty slice
/home/dpirvu/.local/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1671: RuntimeWarning: Degrees of freedom <= 0 for slice.
  keepdims=keepdims)
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:36: MatplotlibDeprecationWarning: savefig() got unexpected keyword argument "rasterize" which is no longer supported as of 3.3 and will become an error in 3.6
In [16]:
if True:
    titls = [r'$\rm Sub-critical \; Bounce$', r'$\rm Sphaleron \; Precursor$']

    home_dor = '/home/dpirvu/project/velocityCOM/paper_thermal/'
    sphaleron = home_dor+'instanton_x1024_phi01.3963_lambda1.5000_fields.dat'
    precursor = home_dor+'precursor_x1024_phi01.3963_lambda1.5000_T0.0900_sim0_fields.dat'

    instanton_sim = extract_data(nLat, sphaleron)
    instanton_sim = instanton_sim[0, :nLat*2//3, nLat//3:2*nLat//3+1][::-1]
    nT, nN = np.shape(instanton_sim)
    print(np.shape(instanton_sim))

    precursor_sim = extract_data(nLat, precursor)
    precursor_sim = precursor_sim[0, :nT, nLat//3:2*nLat//3+1][::-1]
    print(np.shape(precursor_sim))

    exts = np.array([-nN//2,nN//2,0,nT])*dx*np.sqrt(m2(lamb))

    fig, ax = plt.subplots(1, 2, figsize = (8,3))
    im0 = ax[0].imshow(instanton_sim, aspect='auto', interpolation='none', extent=exts, origin='lower', cmap='RdBu')
    clb0 = plt.colorbar(im0, ax = ax[0], ticks = mticker.MultipleLocator(np.pi/2), format = mticker.FuncFormatter(multiple_formatter()))
    clb0.ax.set_title(r'$\bar{\phi}$')

    im1 = ax[1].imshow(precursor_sim, aspect='auto', interpolation='none', extent=exts, origin='lower', cmap='RdBu')
    clb1 = plt.colorbar(im1, ax = ax[1], ticks = mticker.MultipleLocator(np.pi/2), format = mticker.FuncFormatter(multiple_formatter()))
    clb1.ax.set_title(r'$\bar{\phi}$')

    for ai, aa in enumerate(ax):
        aa.text([7.,6.4][ai], [12.,18.][ai], titls[ai], ha='center', va='center', \
                bbox={'boxstyle':'round','facecolor':'white','alpha':0.85,'edgecolor':'k','pad':0.3}, fontsize=10)

        aa.grid(ls=':', color='darkgray', alpha=0.7)
        aa.xaxis.set_minor_locator(MultipleLocator(1))
        aa.yaxis.set_minor_locator(MultipleLocator(1))
        aa.set_xlabel(r'$m \, r$')
        aa.set_ylabel(r'$m \, t$')
        aa.yaxis.set_ticks_position('both')
        aa.xaxis.set_ticks_position('both')
        aa.tick_params(which='both', axis="y", direction="in")
        aa.tick_params(which='both', axis="x", direction="in")
        # a = aa.get_xticks().tolist()[1:-1:]
       # a = [round(al,2) for al in a]
       # aa.set_xticks(a)
        # aa.set_xticklabels(a)
        a = aa.get_yticks().tolist()[1:-1:2]
        a = [round(al,2) for al in a]
        aa.set_yticks(a)
        a = [r'${:.0f}$'.format(al) for al in a]
        aa.set_yticklabels(a)

    plt.savefig('./plots/standard_oscillons.pdf', dpi=500, rasterize=True)
    plt.show()
(682, 342)
(682, 342)
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:51: MatplotlibDeprecationWarning: savefig() got unexpected keyword argument "rasterize" which is no longer supported as of 3.3 and will become an error in 3.6
In [17]:
instanton_sim = instanton_sim - np.mean(instanton_sim)
nT, nN = np.shape(instanton_sim)
aHD_inst = get_HT(instanton_sim)

precursor_sim = precursor_sim - np.mean(precursor_sim)
aHD_prec = get_HT(precursor_sim)

fig, ax = plt.subplots(2, 2, figsize=(8,5))
im1 = ax[0,0].imshow(aHD_inst, aspect='auto', interpolation='none', origin='lower')
im2 = ax[1,0].imshow(aHD_prec, aspect='auto', interpolation='none', origin='lower')

l1 = get_osc_trajectory(aHD_inst, extent=deltae)
l2 = get_osc_trajectory(aHD_prec, extent=deltae)

ylist = np.arange(nT)
ax[0,0].plot(l1, ylist, color='orange', linewidth=1)
ax[1,0].plot(l2, ylist, color='orange', linewidth=1)

avOscillon_inst = instanton_sim[np.arange(nT), l1][::-1]
ax[0,1].plot(np.arange(nT), avOscillon_inst)
avOscillon_prec = precursor_sim[np.arange(nT), l2][::-1]
ax[1,1].plot(np.arange(nT), avOscillon_prec)
ax[1,1].plot(np.arange(nT), avOscillon_inst)
ax[1,1].plot(np.arange(len(avOscillon)), avOscillon)
ax[1,1].plot(np.arange(len(avOS)), avOS)

ax[0,1].set_xlabel(r'$t$')
ax[1,1].set_xlabel(r'$t$')
ax[0,1].set_ylabel(r'$\bar{\phi}(osc,t)$')
ax[1,1].set_ylabel(r'$\bar{\phi}(osc,t)$')
ax[0,1].axhline(0, ls=':', color='k')
ax[1,1].axhline(0, ls=':', color='k')

plt.colorbar(im1, ax = ax[0,0])
plt.colorbar(im2, ax = ax[1,0])
plt.tight_layout()
plt.show()
In [18]:
fig, ax = plt.subplots(1, 1, figsize=(6, 3))
simLists, oscillonAmp, means = oscillonData[:,0], oscillonData[:,1], oscillonData[:,2]
for si, sim in enumerate(simLists):
    oscillon = oscillonAmp[si]
    mn = means[si]
    freqs, PS_oscillon, ind = oscillonFreqData[si]

    nT = len(oscillon)
    plt.loglog(freqs[1:], PS_oscillon[1:], linewidth=0.1, color='k')

y_fft  = np.fft.rfft(avOscillon_inst)
freqs1 = np.fft.rfftfreq(len(avOscillon_inst), d=1.) * len(avOscillon_inst)
psyfd  = (np.conj(y_fft)*y_fft).real
freqmax1 = freqs1[np.argmax(psyfd)]
plt.loglog(freqs1, psyfd, ls='-', marker='o', ms=1, label=r'$\rm Sub-Critical \; Sphaleron$')

y_fft  = np.fft.rfft(avOscillon_prec)
freqs2 = np.fft.rfftfreq(len(avOscillon_prec), d=1.) * len(avOscillon_prec)
psyfd  = (np.conj(y_fft)*y_fft).real
freqmax2 = freqs2[np.argmax(psyfd)]
plt.loglog(freqs2, psyfd, ls='-', marker='o', ms=1, label=r'$\rm Precursor$')

if False:
    avOS, _ = tolerant_mean(oscillonAmp)
    y_fft  = np.fft.rfft(avOS)
    freqs3 = np.fft.rfftfreq(len(avOS), d=1.) * len(avOS)
    psyfd  = (np.conj(y_fft)*y_fft).real
    freqmax3 = freqs3[np.argmax(psyfd)]
    plt.loglog(freqs3, psyfd, ls='-', marker='o', ms=1, label=r'$\rm Average \; Oscillon$')

plt.xlim((freqs2[1], freqs2[-100]))
leg = ax.legend(ncol=1, loc='best', fontsize=10, bbox_to_anchor=(1, 0.8))
plt.grid(True, ls=':', alpha=0.5)

plt.tight_layout()
plt.savefig('./plots/all_freqs.pdf', dpi=500, rasterize=True)
plt.show()
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:36: MatplotlibDeprecationWarning: savefig() got unexpected keyword argument "rasterize" which is no longer supported as of 3.3 and will become an error in 3.6
In [19]:
fig, ax = plt.subplots(1, 1, figsize=(3, 2.5))

plt.axvline(dk/np.sqrt(m2(lamb))*freqmax1, color='r', ls='-.', label=r'$\rm Sub-Critical \; Sphaleron$')
plt.axvline(dk/np.sqrt(m2(lamb))*freqmax2, color='k', ls=':', label=r'$\rm Precursor$')
print(freqmax1, freqmax2)

all_freqs, all_PS_oscillon, all_ind = oscillonFreqData[:,0], oscillonFreqData[:,1], oscillonFreqData[:,2]

freqs_vec = np.zeros(len(all_freqs))
for ii, ind in enumerate(all_ind):
    freqs_vec[ii] = all_freqs[ii][ind]

ah = plt.hist(dk/np.sqrt(m2(lamb)) * np.array(all_ind), bins=50, density=False, label=r'$\rm Oscillons$')
#plt.axvline(dk*freqmax3, color='g', ls=':', label=r'$\rm Average \; Oscillon$')

avfq, error1 = tolerant_mean(oscillonFreqData[:,0])
[plt.axvline(dk/np.sqrt(m2(lamb))*ii, color='k', ls='-', alpha=0.6, linewidth=0.1) for ii in avfq]

avfq *= dk/np.sqrt(m2(lamb))
ax.set_xlim((avfq[1], avfq[10]))

ax.set_ylabel(r'$\rm PDF$')
ax.set_xlabel(r'$\omega_{\rm osc}$')
leg = ax.legend(ncol=1, loc='best', fontsize=10, bbox_to_anchor=(1, 0.8))

plt.savefig('./plots/characteristic_osc_frequency.pdf', dpi=500, rasterize=True)
plt.plot()
5.0 5.0
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:128: RuntimeWarning: Mean of empty slice
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:26: MatplotlibDeprecationWarning: savefig() got unexpected keyword argument "rasterize" which is no longer supported as of 3.3 and will become an error in 3.6
Out[19]:
[]
In [ ]:
 
In [ ]: